Heavy nuclei at the end of the cosmic ray spectrum? 
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We provide an account of the possible acceleration of iron nuclei up to energies ~ 300 EeV in the 
nearby, metally-rich starburst galaxy NGC 253. It is suggested that particles can escape from the 
nuclear region with energies of ~ 10 eV and then could be reaccelerated at the terminal shock 
of the galactic superwind generated by the starburst, avoiding in this way the photodisintegration 
expected if the nuclei were accelerated in the central region of high photon density. We have also 
made estimates of the expected arrival spectrum, which displays a strong dependency with the 
energy cutoff at the source. 

PACS number(s): 96.40 - 98.70.S - 95.85.R - 13.85.T 



I. INTRODUCTION 

The discovery of extensive air shower events with en- 
ergies above 100 EeV (see Ref. |5J] for a recent survey) 
confirms that the cosmic ray (CR) spectrum does not end 
with the expected Greisen-Zatsepin-Kuz'min (GZK) cut- 
off Q . The origin of this GZK cutoff is energy degrada- 
tion of the CR particles (usually assumed to be nucleons 
and nuclei) by resonant scattering processes with the dif- 
fuse background radiation that permeates the universe. 
The observed tail of the spectrum could be, consequently, 
originated in a bunch of nearby sources Q. 

Preferred sites for proton acceleration are astrophys- 
ical scenarios where large-scale shocks occur, as for in- 
stance in the hot spots of powerful radio-galaxies B . It 
could appear that, because of the high-energy cutoff of 
shock acceleration increases with the charge number of 
the nucleus, heavy ions would be nice candidates for ultra 
high-energy CRs. However, AGNs and radio galaxies are 
widely thought as regions of very low metallicity and, in 
addition, it is well established that above 200 EeV nuclei 
should be photodissociated by the 2.73 K photon back- 
ground in a few Mpc p]. Thereupon, nuclei acceleration 
up to the highest energies within astrophysical environ- 
ments is seldom considered in the literature. 

Nonetheless, there has been a recently renewed interest 
in the propagation of heavy nuclei [^]-|f . This renewal is 
mainly sustained by two facts: i) a medium mass nucleus 
is the particle that provides the best fits of the shower 
development of the highest energy CR event Jl0| , and ii) 
the arrival direction of such event roughly points towards 
the nearby metally-rich galaxy M82 

Despite the aforementioned studies on nuclei propaga- 
tion, it is far from clear whether iron or other heavy nuclei 
can be accelerated up to energies ~ 300 EeV in starburst 
galaxies like M82. One could naively expect that, since 
the size scale of the starburst region is of the order of the 
gyroradius of 300 EeV-(Z = 26) nuclei, strong shocks 
could diffusively accelerate these ions to ultra-high ener- 
gies. But this hope fades away as soon as one notices the 
large photon energy densities (mostly in the far infrared) 



measured in the central regions of these kind of galaxies: 
iron nuclei are photodisintegrated long before they can 
reach the required Lorentz factors. 

We shall argue in this paper that, despite the men- 
tioned problem, iron nuclei can be actually acclerated in 
nearby starburst galaxies up to energies ~ 300 EeV if a 
two-step process is involved. The crucial point is that for 
energies above ~ 10 15 eV, acceleration occurrs in the ter- 
minal shock of the starburst superwind, well outside the 
problematic nuclear region. Since NGC 253 is a south- 
ern object which has been scarcely discussed in relation 
to CRs (for a brief discussion of M82 as CR accelerator 
see jy]), we shall focus on it for our quantitative esti- 
mates. However, we emphasize that due to the similarity 
between both galaxies, our conclusions will be correct to 
M82 within the order of magnitude. Let us start with 
a recapitulation of some observational features of NGC 
253. 



II. THE STARBURST GALAXY NGC 253 

Starbursts are galaxies undergoing a massive and large- 
scale star formation episode. Their characteristic signa- 
tures are strong infrared emission (originated in the high 
levels of interstellar extinction), a very strong HII- region- 
type emission-line spectrum (due to a large number of 
O and B-type stars), and a considerable radio emission 
produced by recent supernova remnants. Typically, the 
starburst region is confined to the central few hundreds 
of parsecs of the galaxy, a region that can be easily 10 
or more times brighter than the center of normal spiral 
galaxies. 

NGC 253 has been described as the archetypal star- 
burst galaxy by Rieke et al. |Q , and as a prototype of su- 
perwind galaxy by Heckman and collaborators |l3| . This 
object, whose distance is estimated in the range 2.5 - 3.4 
Mpc |I^,|l5|, has been extensively studied from radio to 
7-ray wavelengths [^6|-^8| . More than 60 idividual com- 
pact radio sources have been detected within the central 
200 pc of the nuclear region of NGC 253 pi most of 
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which are supernova remnants (SNRs) of only a few hun- 
dred years old. According to estimates from observations 
at different frequencies the supernova rate is as high as 
0.2-0.3 yr" 1 §||o|. 

The central ~ 80 pc of the galaxy contain around 
24.000 O and 3.000 red supergiant stars [^(J, in addition 
to the SNRs and numerous HII regions. This means that 
the massive star formation rate is ~ 0.1 M Q yr _1 . Strong 
[Fe II] emission has been also detected with a total [Fe 
II] (1.644 /xm) luminosity of ~ 2.8 x 10 39 erg s _1 , which 
reflects the very rich iron production in the supernovae 
and their associated shocks [ pl[ . 

In the light of such a concentrated activity it is not 
surprising that strong physical, morphological, and kine- 
matic evidence for the existence of a galactic superwind 
has been found in NGC 253 p2| , |T3"| ]. Galactic-scale su- 
perwinds are driven by the collective effect of supernovae 
and massive star winds. The high supernovae rate creates 
a cavity of hot gas (~ 10 8 K) whose cooling time is much 
greater than the expansion time scale. Since the wind 
is sufficiently powerful, it can blow out the interstellar 
medium of the galaxy avoiding to remain trapped as a 
hot bubble. As the cavity expands a strong shock front is 
formed on the contact surface with the cool interestellar 
medium. Shock interactions with low and high density 
clouds produce X-ray continuum and optical line emis- 
sion, respectively, that has been directly observed 0]. In 
addition, kiloparsec regions well outside the disk present 
double emission-line profiles with line splitting of 200-600 
km s^ 1 , a clear evidence of mass motion. The morphol- 
ogy of the optical emission line nebulae indicates that the 
outflowing gas is located along the walls of a cone that 
is limb-brightened, typical of a superwind in a blowout 
phase. 

The shock velocity can reach several thousands of kilo- 
meters per second and ions like iron nuclei can be then 
efficiently accelerated in this scenario up to high energies 
(~ 10 20 eV) by Fermi mechanism as we shall discuss in 
the next section. 



III. CR ACCELERATION AT NGC 253 

We suggest that the iron nuclei acceleration in NGC 
253 occurs through a two-step process. In a first stage, 
ions are diffusively accelerated at single supernova shock 
waves within the nuclear region of the galaxy. Energies 
up to ~ 10 15 eV can be achieved in this step p3| . Fe- 
nuclei are not photodissociated in the process despite the 
star burst's central photon density is much larger than 
that of the Milky Way. The continuum spectrum of NGC 
peaks in the far infrared at ~ 100/im, with a luminosity 
of ~ 3 x 10 10 Lq pj] and a photon energy density of 
Uph ~ 200 eV cm~ 3 ||13[|. For such values is straightfor- 
ward that the nuclei interactions with the blue-shifted 
ambient photons are quite below the photodisintegration 
threshold. Energy losses through pair production are also 



negligible. 

By other hand, interactions with the rich interstellar 
gas in the center of the galaxy might disintegrate the ions 
if the escape from the starburst region is dominated by 
diffussion. Multiwavelength molecular line observations 
of the central region show that the average density of the 
molecular clouds is ~ 10 5 cm -3 , with a filling factor < 
10 -3 P5[ . This means an average gas density (mainly H2) 
in the active region within the range 30 - 300 cm -3 Jl7| . 
The cross section for Fe— H2 interactions [a 1463.7 
mb) can be estimated from Rudstam's parametrization 
of proton induced spallation, a = 50 A^ 3 mb (where A t is 
the mass number of the target nucleus) |2(| . The mean 
free path for the iron nuclei is then A ~ (wcpc-Ha) -1 
which results in the range 738 - 7380 pc when the upper 
and lower limits of particle density are considered. The 
central starburst region can be modeled as a disk of 70 pc 
thick with a radius R w 300 pc jl7| . Since the gyroradius 
of the Fe nuclei is ~ 10 -3 pc, they certainly cannot be 
driven out from the starburst by diffusion. 

However, due to the nature of the central region in 
NGC 253, the escape of the iron nuclei is expected to be 
dominated by convection. In fact, the presence of several 
tens of young SNRs with very high expansion velocities 
(~ 12000 km s _1 |p7| ) and thousands of massive O stars 
(with stellar winds of terminal velocities up to 3000 km 
s _1 ) must generate collective plasma motions of several 
thousands of km per second. Then, due to the coupling 
of the magnetic field to the hot plasma, the magnetic 
field is also lifted outwards and forces the cosmic ray gas 
to stream along from the starburst region. 

The relative importance of convection and diffusion in 
the escape of the cosmic rays from a region of disk scale 
height h is given by the dimcnsionlcss parameter, 



where Vq is the convection velocity and n is the cosmic 
ray diffusion coefficient inside the starburst p8| |. When 
q < 1, the cosmic ray outflow is "difussion dominated", 
whereas when q > 1 it is "convection dominated" . As- 
suming for the central region of NGC 253 a convection 
velocity of the order of the expanding SNR shells (i.e. ~ 
10000 km s _1 , a scale height h ~ 35 pc, and a reasonable 
value for the diffusion coefficient n$ ~ 5 x 10 26 cm 2 s" 1 
p9| , we get q ~ 216 and convection dominates the escape 
of the particles. The residence time of the iron nuclei in 
the starburst results ^res ~ h/Vo ~ 1 x 10 11 s. Most of 
the nuclei then escape through the disk in opposite direc- 
tions along the symmetry axis of the system, being the 
total path traveled substantially shorter than the mean 
free path. 

Once the nuclei escape from the central region of the 
galaxy with energies of ~ 10 15 eV, they are injected into 
the galactic-scale wind and experience further accelera- 
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tion at its terminal shockj^] The scale length of this shock 
is of the order of several tens of kpc (see Ref. [H ), so it 
can be considered as locally plane for calculations. The 
shock velocity v s h can be estimated from the empirically 
determined superwind kinetic energy flux E sw and the 
mass flux M generated by the starburst through: 



En- 



(2) 



The shock radius can be approximated by r « v s hT, 
where r is the starburst age. Since the age is about a 
few tens of million years, the maximum energy attain- 
able in this configuration is constrained by the limited 
acceleration time arisen from the finite shock's lifetime. 
For this second step in the acceleration process, the pho- 
ton field energy density drops to values of the order of 
the cosmic background radiation (we are now far from 
the starburst region), and consequently, as we shall dis- 
cuss in the next section, iron nuclei are again safe from 
photodissociation while energy increases from ~ 10 15 to 
10 20 eV. 

In order to estimate the maximum energy that can be 
reached by the nuclei in the second stage of the accelera- 
tion process let us consider the superwind terminal shock 
propagating in a homogeneous medium with an average 
magnetic field B. If we work in the frame where the 
shock is at rest, the upstream flow velocity will be Vi 
(|vi| = u s h) and the downstream velocity, v 2 . The mag- 
netic field turbulence is assumed to lead to isotropiza- 
tion and consequent diffusion of energetic particles which 
then propagate according to the standard transport the- 
ory (3^] . The acceleration time scale is then [^2) : 
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where k is the upstream diffusion coefficient which can 
be written in terms of perpendicular and parallel compo- 
nents to the magnetic field, and the angle 9 between the 
(upstream) magnetic field and the direction of the shock 
propagation: 



Ku cos 



k± sm 



(4) 



Since strong turbulence is expected from the shock we 
can take the Bohm limit for the upstream diffusion coef- 
ficient parallel to the field, i.e. 
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where B\ is the strength of the pre-shock magnetic field 
and E is the energy of the Z-ion. For the k± component 



*CR acceleration at superwind shocks was firstly proposed 
by Jokipii and Morfill during the 80s J3(J in the context of 
our own Galaxy. 



we shall assume, following Biermann |33|, that the free 
mean path perpendicular to the magnetic field is inde- 
pendent of the energy and has the scale of the thickness 
of the shocked layer (r/3). Then, 



-r(vi - v 2 ) 



or, in the strong shock limit, 
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Since the upstream time scale is t a 
rewrite Eq. (H) as: 
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r/(3ui), we 
(8) 



and then, using r — v±t and transforming to the ob- 
server's frame, we get: 
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-ZeBv sh r, 



Bmax « -ZeB— t-T, 
2 M 



(9) 
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in terms of parameters that can be determined from ob- 
servations. 

The predicted kinetic energy and mass fluxes of the 
starburst of NGC 253 derived from the measured IR lu- 
minosity are 2 x 10 42 erg s _1 and 1.2 M Q yr _1 , respec- 
tively p3[ . The starburst age is estimated from numer- 
ical models that use theoretical evolutionary tracks for 
individual stars and make sums over the entire stellar 
population at each time in order to produce the galaxy 
luminosity as a function of time jjij. Fitting the ob- 
servational data these models provide a range of suitable 
ages for the starburst phase that, in the case of NGC 253, 
goes from 5 x 10 7 to 1.6 x 10 s yr (also valid for M82) Q. 
These models must assume a given initial mass function 
(IMF), which usually is taken to be a power-law with 
a variety of slopes. Recent studies has shown that the 
same IMF can account for the properties of both NGC 
253 and M82 ||l[]. Here we shall assume a conservative 
age r = 50 Myr. Finally, the radio and 7-ray emission 
from NGC 253 are well matched by models with B ~ 50/x 
G p^t . With these figures we obtain a maximum energy 
for iron nuclei of: 



E, 



Fc 



3.4 x 10 2 



eV, 



(11) 



a value quite lower (more than two orders of magnitude) 
than the limit imposed by the synchrotron losses |3q] . 
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1. Effective energy loss time for iron nuclei photodis- 



integration through the 2.73K and IR photons. Adapted from 
Epele and Roulet Ref . [7] . 

Regarding the spectral slope of the particles, it is ex- 
pected that they emerge from the central region of NGC 
253, where there are strong stellar winds from the mas- 
sive star population, with an energy index 7 = 2.4 [ |33| , 
which is in fine accordance with the observed radio spec- 
trum oc v~ - 7 observed in the individual SNRs within the 
starburst . After the particles left the nuclear region, 
diffusive losses tend to steepen the spectrum at the termi- 
nal shock of the galactic superwind, where reacceleration 
takes place. The final index will depend on some not well 
known parameters as the Mach number of the terminal 
shock. In what follows, we shall study the propagation 
of the nuclei from NGC 253 to the Earth in models with 
7 = 2.4, 2.5, and 2.6. 



IV. PROPAGATION EFFECTS 

The basic interactions between the universal back- 
ground radiation and nuclei of extremely high energy 
have been first discussed in detail by Puget et al [ |36| . 
Heavy nuclei with energies above a few EeV get atten- 
uated mainly by photodisintegration off the microwave 
background radiation (MBR) and the intergalactic in- 
frared (IR) background photons .0 

The photodisintegration process is dominated by a 
broad maximum designated as the giant dipole resonance 
which peaks in the 7-ray energy range of 10 to 30 MeV 
(nucleus rest frame). In the initial absorption process, 
the photon energy may be given to a single nucleon (al- 
though almost always is shared with others), so the decay 



t Actually, the pair creation process due to interactions with 
the MBR as well as disintegration with the optical background 
photons also attenuate the nuclei energy. However, these pro- 
cesses are not essential for the discussion presented in this 
paper. For details of these interactions the reader is referred 
to Ref. H. 



process might involve the emission of one (or more) nu- 
cleons. The disintegration rate of a nucleus of mass A 
with the subsequent production of i nucleons reads [B7[, 



dw 



n(w) 



2Tw 



dw r w r OAi(w r ) (12) 



where n(w) is the density of photons with energy w in 
the system of reference in which the MBR is at 2.73 K. In 
the formula, w r is the energy of the photons in the rest 
frame of the nucleus, F is the Lorentz factor, and oai 
is the cross section for the interaction. The total cross 
section of photonuclear interaction is well known [p8[ . 
As usual, the MBR is described by Plankian spectrum 
of 2.73 K. For the diffuse IR radio background (photon 
energies between 2 x 10 -3 - 8 x 10 _1 eV) we shall use the 
estimate given in Ref. [p9], 



— = 1.1 x 10~ 4 ( — 
dw VeV/ 



cm _3 eV _1 . 



(13) 



With these figures at hand it is straightforward to com- 
pute the energy loss time of iron nuclei (displayed in 
Fig. 1). It is clear from these calculations that the 
universal background radiation would not affect iron ac- 
celeration in NGC 253 up to a few thousands of EeV. 
Moreover, since the nearness of the starburst galaxy, just 
the interaction with the MBR becomes relevant in their 
travel to Earth. The fractional energy loss as a function 
of the MBR energy (Lorentz factor) have been already 
parametrized 0, 



R(T) = 3.25 x 10- b r-°- b4ci exp(-2.15 x 10 lo /r) s" 



if T e [1. x 10 9 ,3.68 x 10 10 ], and 

i?(r) = i.59xio- 12 r-°- 0698 



(14a) 



(14b) 



if r e [3.68 x 10 10 ,1. x 10 u ]. It is noteworthy that 
t he p ossibl e dis integration histories computed using Eqs. 
( |l4a| ) and (|l4b|) are in very good agreement with Monte 
Carlo simulations when the sources are located near the 
Earth (distances < 10 Mpc) @. 

Using the formalism sketched in ||, it is easily 
obtained the evolution of the differential spectrum 
Q(E g ,t) = KEg 7 S(t — to) of the iron nuclei injected by 
NGC 253 at to. The number of surviving fragments with 
energy E at time t is given by, 



KE~"i +1 
N(E,t)dE = | dE, 



(15) 



where E g denotes the energy at which the nuclei were 
emitted from the source, related to the energy detected 



on Earth by E — E„ e 



-R(T) t/56 



(recall that the Lorentz 
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factor, r = E, 
propagation) 



56, does not result modified during the 
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FIG. 2. Relation between the injection and arrival energies 
for possible propagation distances from NGC 253 according 
to current data. 

In Fig. 2 it is shown the energy degradation for dif- 
ferent flying timescales arising from different propagation 
regimes. The energy spectrum of the surviving fragments 
is degenerate. For instance, for a propagation distance of 
3 Mpc, the composition of the arrival nuclei changes from 
A = 55 (for T w 3 x 10 9 ) to A = 44 (for r w 6 x 10 9 ). If a 
maximum energy of 340 EeV is attainable in the source, 
the relation between the injection and arrival energies is 
a monotonously decreasing function. At this stage, we 
conveniently introduce the modification factor 77, defined 
as the ratio between the modified spectrum and the un- 
modified one. Notice that once the nuclei have energy 
enough to undergo photodisintegration through the gi- 
ant dipole resonance, the value of the modification factor 
is always less than one. 

Contrariwise, if the injection energy has an upper cut- 
off at E > 560 EeV (this can be achieved, for instance, 
with a higher value of r), the function which relates the 
injection and arrival energies becomes multivalued, yield- 
ing a bump- like feature in the modification factor. To un- 
derstand this behavior, recall that nuclei suffer a violent 
disintegration at these energies (see Fig. 1), in such a way 
that, for a Lorentz factor r « 8 x 10 9 , the composition of 
the arrival nuclei drops to A m 33. So, particles injected 



with different energies might arrive with the same energy, 
piling up around 250 EeV, just before the expected cut- 
off. The whole effect can be clearly appreciated in Fig. 
3, where we have plotted the modification factor for two 
different cutoffs in the injection spectrum. 

For longer propagation distances, the pile-up shifts to 
lower energies (e.g. to « 200 EeV for 3.4 Mpc), broaden- 
ing its profile. On the other hand, for shorter distances 
to the nuclei-emitting source there is no spectral bump, 
i.e., the spectrum is monotonic. Thus, the probability of 
detecting events diminishes with increasing Lorentz fac- 
tor. 

We remark that the change of the spectral index does 
not significantly affect the shape of the modification fac- 
tor. The height of the pile-up diminishes with the steep- 
ening of the injected spectrum, but this effect is accom- 
panied by a simultaneous drop in the CR flux of the pre- 
ceeding bin of energy, resulting in the same overall shape 
for rj, although downshifted (see Figs. 4 and 5 of ref M). 




FIG. 3. Modification factor for initial iron nuclei from NGC 
253, assuming a differential power law injection spectrum with 
spectral index 7 = 2.4. Solid (dotted) line stands for the case 
with an upper cutoff at the injection energy of 560 EeV (340 
EeV). 

Events from the pile-up are about 50% more proba- 
ble to be detected than those at energies immediately 
lower. The reason is that nuclei with energies > 560 EeV 
would be almost completely photodisintegrated during 
their journey to Earth, in such a way that the surviving 
fragments end at energies of the pile-up, changing the 
relative detection probabilities. Thus, we have the inter- 
esting result that the existence of a sufficient high cutoff 
(at Lorentz factors, say, T > 8 x 10 9 ) in the source accel- 
eration mechanism favors the detection of events around 
250 EeV from the starburst galaxy NGC 253. 



*Notice that the (Z, ^-dependence of R is roughly cancelled 
by dividing by A in the exponent. This implies that R can in 
fact be integrated down to lower, spallated A values and still 
be reasonably accurate |p6[. 



V. CONCLUDING REMARK 

Might heavy nuclei be primaries at the end of the cos- 
mic ray spectrum? If the spectrum extends over the up- 
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to-now observed energies, the answer will be certainly 
"no" . Heavy nuclei with energies above 200 EeV could 
not propagate for more than 10 Mpc. Besides, Lorentz 
factors above 2 x 10 10 require large astrophysical regions 
for acceleration where the relic photon density would be 
sufficient to provide a significant loss mechanism for the 
nuclei (see Fig. 1). 

However, if the highest cosmic ray energy is of the or- 
der of 300 EeV, heavy nuclei accelerated in the nearby 
starburst galaxies (NGC 253 and M82) can be present at 
the end of the spectrum. We have shown that these nu- 
clei, originated in the central regions of the galaxies, can 
be accelerated without suffering catastrophic interactions 
in a two-step process that involves supernova remnant 
shock waves, and the large scale terminal shock produced 
by the superwind that flows from the starbursts. 

Anchordoqui et al. J8| have suggested that the lack of 
data at energies immediately lower than the two "super- 
GZK events" recorded to date, could be the result of a 
different primary composition of disintegrated iron nuclei 
at the end of the spectrum. This speculation might find 
some support from our calculations since they suggest 
that the production of the energetic nuclei in next-door 
galaxies (NGC 253 and M82) is at least feasible. In any 
case, we have made some concrete predictions that could 
be tested with the forthcoming facilities of the South- 
ern Auger Observatory M] , and the satellite experiment 
called OWL 0. 
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